LAB10 data
RStudio
CloudCompare
ArcGIS Pro
So we’ve made it this far into the quarter and it’s time for the last lab. There is not final exam but this lab will consist of many of the past lessons from the previous 9. By now you should have a repository of code and lidar processing examples in the form of your past assignments. Most of this final lab comes from those past assignments and reuses their code. So there will be more limited examples here meaning less copy/paste.
The data in this lab also consists of past lab examples except for
some extra TLS data from Pack Forest and the Hall of Mosses in the Hoh
Rainforest: 2019-08-18 Pack Forest TMLS.laz&
2019-08-22 Hall of Mosses TMLS.laz. You will need to
download extra data from the lidar portal in Part 1.
You will have 2 plots for this part centered at:
Pack Forest - 1185950 555180 (EPSG 2927) / -122.31625, 46.84141 (WGS 84):
ONP - 797785 941189 (EPSG 2927) / -123.93314, 47.86372 (WGS 84):
Download the ALS point clouds at those two locations from the Washington DNR lidar portal website.
Read in the two ALS laz files and the two TLS laz files into R
ALS files (you downloaded from the DNR) (Assuming Pack 2013)
TLS files (provided via canvas):
PFTLS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/2019-08-18 Pack Forest TMLS.laz')
HMTLS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/2019-08-22 Hall of Mosses TMLS.laz')
PFALS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/WA_031_rp.laz')
HMALS <- readLAS('Lab 10/Lab_10_2023/LAB10_Data_2023/q47123g8201.laz')
crs in R, we need them all to be in the same
crsNEWCRS <- crs("EPSG:2927")
crs(PFTLS) <- NEWCRS
crs(HMTLS) <- NEWCRS
crs(PFALS) <- NEWCRS
crs(HMALS) <- NEWCRS
clip_circle the files to the plot extent:
?clip_circle
PFTLSc <- clip_circle(PFTLS,1185950,555180,50)
HMTLSc <- clip_circle(HMTLS,797785,941189,50)
PFALSc <- clip_circle(PFALS,1185950,555180,50)
HMALSc <- clip_circle(HMALS,797785,941189,50)
Future steps for the TLS will only use these clips
writeLAS the four lidar clips and bring them into
CloudCompare to produce a figure.
?writeLAS
writeLAS(PFTLSc)
writeLAS(HMTLSc)
writeLAS(PFALSc)
writeLAS(HMALSc)
REQUIRED FIGURES: An image of the ALS and TLS plotted together for both the Pack Forest (PF) and Hall of Mosses (HM) sites (two screen shots). Example given has the TLS cloud in gray scale and the ALS cloud in multicolor for the Hall of Mosses forest plot.
Now back in R:
pmf or csf and the
classify_ground command.mycsf <- csf(FALSE, 1, 1, rigidness = 3)
#OR
mypmf <- pmf(ws = 0.5, th = 0.1)
PFTLSg <- classify_ground(PFTLSc, mycsf)
HMTLSg <- classify_ground(HMTLSc, mycsf)
plot(PFTLSg, color = "Classification")
rasterize_terrain to make DTMs for both the TLS
point cloud clipsPFDTM1 <- rasterize_terrain(PFTLSg, res = 1, algorithm = knnidw())
HMDTM1 <-grid_terrain(HMTLSg, res = 1, algorithm = knnidw())
## Warning: There were 1 degenerated ground points. Some X Y Z coordinates were
## repeated. They were removed.
## Warning: There were 172 degenerated ground points. Some X Y coordinates were
## repeated but with different Z coordinates. min Z were retained.
plot_dtm3d to make sure your DTMs are good quality. You
will likely need to try different values in the production of the
DTMs.plot_dtm3d(PFDTM1)
plot_dtm3d(HMDTM1)
REQUIRED FIGURES: Two screen shots of your two produced dtms. Fully captioned with the technique used to make them.
Normalize your point clouds using the DTMs from the TLS.
normalize_height()PFTLScN <- normalize_height(PFTLSc,PFDTM1)
HMTLScN <- normalize_height(HMTLSc,HMDTM1)
PFALScN <- normalize_height(PFALSc,PFDTM1)
HMALScN <- normalize_height(HMALSc,HMDTM1)
Create a DSM for all four clips using
rasterize_canopy, Your choice what algorithm to use.
?rasterize_canopy
PFTLS_DSM <- rasterize_canopy(PFTLScN,res = 1.5, p2r(2))
PFALS_DSM <- rasterize_canopy(PFALScN,res = 1.5, p2r(2))
HMTLS_DSM <- rasterize_canopy(HMTLScN,res = 1.5, p2r(2))
HMALS_DSM <- rasterize_canopy(HMALScN,res = 1.5, p2r(2))
plot_dtm3d(PFTLS_DSM)
plot(PFALScN)
plot(HMTLScN)
plot(HMALScN)
?rumple_index
rumple_index(PFTLS_DSM)
## [1] 4.929346
rumple_index(PFALS_DSM)
## [1] 5.576009
rumple_index(HMTLS_DSM)
## [1] 3.509282
rumple_index(HMALS_DSM)
## [1] 5.997609
REQUIRED FIGURES: Four screen shots. One of each of your produced DSMs with the rumple index included in the caption along with the algorithm used to produce the DSM.
- Use
focal statistics to smooth
your DSM. You can do this step before generating the screenshots and
rumple index but the smoothing may cause an error with the rumple
index.
?focal
## Help on topic 'focal' was found in the following packages:
##
## Package Library
## raster /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
## terra /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
##
##
## Using the first match ...
PFTLS_DSMs <- focal(PFTLS_DSM, w = matrix(1,3,3), fun = mean)
PFALS_DSMs <- focal(PFALS_DSM, w = matrix(1,3,3), fun = mean)
HMTLS_DSMs <- focal(HMTLS_DSM, w = matrix(1,3,3), fun = mean)
HMALS_DSMs <- focal(HMALS_DSM, w = matrix(1,3,3), fun = mean)
plot_dtm3d(PFTLS_DSMs)
plot_dtm3d(PFALS_DSMs)
plot_dtm3d(HMTLS_DSMs)
plot_dtm3d(HMALS_DSMs)
segment_trees using the normalized las clips and the
smoothed DSMs.PFTLStrees <- segment_trees(PFTLScN, lidR::watershed(PFTLS_DSMs, th = 10))
plot(PFTLStrees,color = "treeID", pal = pastel.colors)
PFALStrees <- segment_trees(PFALScN, lidR::watershed(PFALS_DSMs, th = 10))
plot(PFALStrees,color = "treeID", pal = pastel.colors)
HMTLStrees <- segment_trees(HMTLScN, lidR::watershed(HMTLS_DSMs, th = 10))
plot(HMTLStrees,color = "treeID", pal = pastel.colors)
HMALStrees <- segment_trees(HMALScN, lidR::watershed(HMALS_DSMs, th = 10))
plot(HMALStrees,color = "treeID", pal = pastel.colors)
#OR
PFTLSttops <- locate_trees(PFTLS_DSMs, lmf(ws = 3.28, hmin = 6.56))
algo1 <- dalponte2016(PFTLS_DSMs, PFTLSttops, th_tree = 0.1, th_seed= 0.1, max_cr = 10, th_cr = 0.1)
algo2 <- silva2016(PFTLS_DSMs, PFTLSttops, max_cr_factor = 0.6, exclusion = 0.3)
PFTLStrees1 <- segment_trees(PFTLScN, algo1)
plot(PFTLStrees1,color = "treeID", pal = pastel.colors)
PFTLStrees2 <- segment_trees(PFTLScN, algo2)
plot(PFTLStrees2,color = "treeID", pal = pastel.colors)
PFALStrees2 <- segment_trees(PFALScN, lidR::watershed(PFALS_DSMs, th = 10))
plot(PFALStrees,color = "treeID", pal = pastel.colors)
HMTLStrees <- segment_trees(HMTLScN, lidR::watershed(HMTLS_DSMs, th = 10))
plot(HMTLStrees,color = "treeID", pal = pastel.colors)
HMALStrees <- segment_trees(HMALScN, lidR::watershed(HMALS_DSMs, th = 10))
plot(HMALStrees,color = "treeID", pal = pastel.colors)
REQUIRED FIGURES: Four screen shots: for the ALS & TLS in each plot of the lidar data colored by the tree segmentation for each of the normalized lidar clips. Fully captioned with the algorithms used to create the tree segmentation.
Calculate cloud metrics on the normalized point cloud for the ALS and TLS above 2m height (6.56ft)
?cloud_metrics
PFTLS_cm <- cloud_metrics((filter_poi(PFTLScN,Z > 6.56)), ~stdmetrics_z(Z))
PFTLS_cm
## $zmax
## [1] 140.6
##
## $zmean
## [1] 25.78389
##
## $zsd
## [1] 28.6307
##
## $zskew
## [1] 1.809775
##
## $zkurt
## [1] 4.91211
##
## $zentropy
## [1] 0.7387645
##
## $pzabovezmean
## [1] 22.80502
##
## $pzabove2
## [1] 100
##
## $zq5
## [1] 7
##
## $zq10
## [1] 7.5
##
## $zq15
## [1] 8
##
## $zq20
## [1] 8.5
##
## $zq25
## [1] 9
##
## $zq30
## [1] 9.7
##
## $zq35
## [1] 10.4
##
## $zq40
## [1] 11.1
##
## $zq45
## [1] 12
##
## $zq50
## [1] 12.9
##
## $zq55
## [1] 13.9
##
## $zq60
## [1] 15.3
##
## $zq65
## [1] 16.9
##
## $zq70
## [1] 19
##
## $zq75
## [1] 23
##
## $zq80
## [1] 32.2
##
## $zq85
## [1] 61.8
##
## $zq90
## [1] 80.6
##
## $zq95
## [1] 97.4
##
## $zpcum1
## [1] 55.57681
##
## $zpcum2
## [1] 78.54518
##
## $zpcum3
## [1] 82.06199
##
## $zpcum4
## [1] 83.66011
##
## $zpcum5
## [1] 86.83228
##
## $zpcum6
## [1] 91.50673
##
## $zpcum7
## [1] 95.48767
##
## $zpcum8
## [1] 99.00548
##
## $zpcum9
## [1] 99.97011
PFALS_cm <- cloud_metrics((filter_poi(PFALScN,Z > 6.56)), ~stdmetrics_z(Z))
PFALS_cm
## $zmax
## [1] 148.57
##
## $zmean
## [1] 90.02078
##
## $zsd
## [1] 42.76528
##
## $zskew
## [1] -0.9650657
##
## $zkurt
## [1] 2.296773
##
## $zentropy
## [1] 0.8806774
##
## $pzabovezmean
## [1] 69.47884
##
## $pzabove2
## [1] 100
##
## $zq5
## [1] 10.04
##
## $zq10
## [1] 15.74
##
## $zq15
## [1] 19.4475
##
## $zq20
## [1] 25.13
##
## $zq25
## [1] 73.0275
##
## $zq30
## [1] 89.14
##
## $zq35
## [1] 96.9475
##
## $zq40
## [1] 102.82
##
## $zq45
## [1] 106.3125
##
## $zq50
## [1] 109.495
##
## $zq55
## [1] 111.9375
##
## $zq60
## [1] 114.04
##
## $zq65
## [1] 116.15
##
## $zq70
## [1] 118.59
##
## $zq75
## [1] 120.63
##
## $zq80
## [1] 122.53
##
## $zq85
## [1] 125.02
##
## $zq90
## [1] 127.975
##
## $zq95
## [1] 131.5575
##
## $zpcum1
## [1] 8.970482
##
## $zpcum2
## [1] 22.23182
##
## $zpcum3
## [1] 23.87329
##
## $zpcum4
## [1] 24.13247
##
## $zpcum5
## [1] 25.19798
##
## $zpcum6
## [1] 30.0216
##
## $zpcum7
## [1] 41.79986
##
## $zpcum8
## [1] 70.59755
##
## $zpcum9
## [1] 96.74586
HMTLS_cm <- cloud_metrics((filter_poi(HMTLScN,Z > 6.56)), ~stdmetrics_z(Z))
HMTLS_cm
## $zmax
## [1] 110.178
##
## $zmean
## [1] 28.15648
##
## $zsd
## [1] 20.19887
##
## $zskew
## [1] 1.026229
##
## $zkurt
## [1] 3.162441
##
## $zentropy
## [1] 0.8600562
##
## $pzabovezmean
## [1] 37.9027
##
## $pzabove2
## [1] 100
##
## $zq5
## [1] 7.292
##
## $zq10
## [1] 8.09
##
## $zq15
## [1] 9.077
##
## $zq20
## [1] 10.4
##
## $zq25
## [1] 11.863
##
## $zq30
## [1] 13.384
##
## $zq35
## [1] 14.856
##
## $zq40
## [1] 16.642
##
## $zq45
## [1] 18.863
##
## $zq50
## [1] 21.21
##
## $zq55
## [1] 23.4359
##
## $zq60
## [1] 26.506
##
## $zq65
## [1] 30.843
##
## $zq70
## [1] 36.105
##
## $zq75
## [1] 40.834
##
## $zq80
## [1] 45.895
##
## $zq85
## [1] 52.703
##
## $zq90
## [1] 60.056
##
## $zq95
## [1] 69.008
##
## $zpcum1
## [1] 22.20074
##
## $zpcum2
## [1] 51.81598
##
## $zpcum3
## [1] 66.94687
##
## $zpcum4
## [1] 78.20104
##
## $zpcum5
## [1] 86.71497
##
## $zpcum6
## [1] 93.61396
##
## $zpcum7
## [1] 97.68797
##
## $zpcum8
## [1] 99.47636
##
## $zpcum9
## [1] 99.97858
HMALS_cm <- cloud_metrics((filter_poi(HMALScN,Z > 6.56)), ~stdmetrics_z(Z))
HMALS_cm
## $zmax
## [1] 111.88
##
## $zmean
## [1] 68.30136
##
## $zsd
## [1] 23.11787
##
## $zskew
## [1] -0.7348024
##
## $zkurt
## [1] 2.969791
##
## $zentropy
## [1] 0.9386929
##
## $pzabovezmean
## [1] 56.63975
##
## $pzabove2
## [1] 100
##
## $zq5
## [1] 20.18
##
## $zq10
## [1] 34.552
##
## $zq15
## [1] 43.418
##
## $zq20
## [1] 49.264
##
## $zq25
## [1] 54.43
##
## $zq30
## [1] 59.81
##
## $zq35
## [1] 63.262
##
## $zq40
## [1] 66.23
##
## $zq45
## [1] 69.464
##
## $zq50
## [1] 72.74
##
## $zq55
## [1] 75.05
##
## $zq60
## [1] 77.552
##
## $zq65
## [1] 80.428
##
## $zq70
## [1] 83.04
##
## $zq75
## [1] 85.66
##
## $zq80
## [1] 88.512
##
## $zq85
## [1] 91.902
##
## $zq90
## [1] 95.058
##
## $zq95
## [1] 98.868
##
## $zpcum1
## [1] 2.131886
##
## $zpcum2
## [1] 5.826271
##
## $zpcum3
## [1] 9.666314
##
## $zpcum4
## [1] 16.19439
##
## $zpcum5
## [1] 26.20498
##
## $zpcum6
## [1] 41.56515
##
## $zpcum7
## [1] 61.28178
##
## $zpcum8
## [1] 81.43538
##
## $zpcum9
## [1] 96.875
PFTLS_cma <- cloud_metrics(PFTLScN, ~sum(Z>6.56)/sum(Z>-1))
PFTLS_cma
## [1] 0.5067569
PFALS_cma <- cloud_metrics(PFALScN, ~sum(Z>6.56)/sum(Z>-1))
PFALS_cma
## [1] 0.963117
HMTLS_cma <- cloud_metrics(HMTLScN, ~sum(Z>6.56)/sum(Z>-1))
HMTLS_cma
## [1] 0.5147704
HMALS_cma <- cloud_metrics(HMALScN, ~sum(Z>6.56)/sum(Z>-1))
HMALS_cma
## [1] 0.6470487
| name | Zmax | Zmean | Zsd | Zq25 | Zq95 | CanopyCover |
|---|---|---|---|---|---|---|
| Pack Forest ALS | 148.570 | 90.02078 | 42.76528 | 73.0275 | 131.5575 | 0.9631170 |
| Pack Forest TLS | 140.600 | 25.78389 | 28.63070 | 9.0000 | 97.4000 | 0.5067569 |
| ONP ALS | 111.880 | 68.30136 | 23.11787 | 54.4300 | 98.8680 | 0.6470487 |
| ONP TLS | 110.178 | 28.15648 | 20.19887 | 11.8630 | 69.0080 | 0.5147704 |
That is all for the processing of the ALS/TLS plot
This part is dealing with some more recent coding so I’ll provide it here. The questions will refer to the outputs
You will be comparing ALS cloudmetrics and Canopy Height Models (CHMs) to GEDI for the Hall of Mosses area of the Olympic National Park
For the purpose of this lab, we are going to assume that there is no error in the GEDI pulse locations. In reality there is some error and that can account for oddities in the data.
To complete this portion, of the lab you will need to:
Download The DTM and DSM from the Olympic Peninsula
Read them into R and assign the ALS the correct crs
ONP_DTM <- rast('Lab 10/Lab_10_2023/LAB10_Data_2023/hoh_2013_dtm_4.tif')
ONP_DSM <- rast ('Lab 10/Lab_10_2023/LAB10_Data_2023/hoh_2013_dsm_4.tif')
Subtract the DTM from the DSM for an easy Canopy Height Model
#ONPCHM <- ONPDSM-ONPDTM
#plot(ONPCHM)
ONPCHM <- rast('Lab 10/Lab_10_2023/LAB10_Data_2023/out/ONPCHM.tif')
#writeRaster(ONPCHM, ".....................ONPCHM.tif")
Now read in GEDI level 2a data from the data folder for this lab
?getLevel2AM
gedilevel2a<-readLevel2A(level2Apath = file.path ("Lab 10/Lab_10_2023/LAB10_Data_2023/HM_level2a_clip.h5"))
Compute a series of statistics of GEDI Level2A RH100 metric
#Get GEDI Elevation and Height Metrics (GEDI Level2A)
?getLevel2AM
level2AM<-getLevel2AM(gedilevel2a)
head(level2AM[,c("beam","shot_number","elev_highestreturn","elev_lowestmode","rh100")])
Convert shot_number as “integer64” to “character” and the make the
GEDI Level2AM points into a spatial object with sf
level2AM$shot_number<-as.character(level2AM$shot_number)
level2AM_sf<-st_as_sf(level2AM, coords = c("lon_lowestmode", "lat_lowestmode"), crs = 4326)
Define a function for computing GEDI stats
mySetOfMetrics = function(x)
{
metrics = list(
min =min(x), # Min of x
max = max(x), # Max of x
mean = mean(x), # Mean of x
sd = sd(x)# Sd of x
)
return(metrics)
}
Create a RH100 raster with 0.005 res in degrees of lat / long.
?gridStatsLevel2AM
rh100metrics<-gridStatsLevel2AM(level2AM = level2AM, func=mySetOfMetrics(rh100), res=0.005)
crs(rh100metrics) <- "EPSG:4326" #same as GEDI
rh100metrics_prj <- terra::project(rast(rh100metrics), "EPSG:2927") #transform to WA state plane south
plot(ext(770000, 850000, 900000, 970000), type = "lines", buffer = TRUE)
plot(rh100metrics_prj$max, axes = F, add = TRUE, legend =T)
plot(ext(770000, 850000, 900000, 970000), type = "lines", buffer = TRUE)
plot(ONPCHM, axes = F, add = TRUE, legend =T)
writeRaster(rh100metrics_prj, "Lab 10/Lab_10_2023/LAB10_Data_2023/out/rh100metrics_prj.tif")
REQUIRED FIGURES: Create 2 figures: a map from the ONPCHM and a map of the same area using the GEDI RH100 max data. Write out the ONPCHM and GEDI RH100 to open in ArcGIS Pro and put them next to each other for comparison and qualitatively assess the spatial patterns. In your caption explain if they generally agree or disagree with each other
Extract ONPCHM elevation values with the new projected GEDI Level2A points
bind=TRUE for combining extracted elevation values
to the points dataframeONPCHM <- rast("Lab 10/Lab_10_2023/LAB10_Data_2023/out/ONPCHM.tif")
level2AGeo_sf_prj <- st_transform(level2AM_sf, "EPSG:2927") #transform to WA state plane south
ONP_CHM_GEDI<- (terra::extract(ONPCHM, level2AGeo_sf_prj, bind = TRUE))
ONP_CHM_GEDI_FILTER <- ONP_CHM_GEDI[-c(is.na(ONP_CHM_GEDI$hoh_2013_dsm_4)),] #removing NA values
ONP_CHM_GEDI_FILTER <- ONP_CHM_GEDI_FILTER[ONP_CHM_GEDI_FILTER$rh100 > 0,] #lets remove the weird GEDI rh100 values that are 0
mapview(st_as_sf(ONP_CHM_GEDI_FILTER))
Create a regression model between GEDI and ALS CHM
reg <- lm(ONP_CHM_GEDI_FILTER$hoh_2013_dsm_4 ~ (ONP_CHM_GEDI_FILTER$rh100))
summary(reg)
plot((ONP_CHM_GEDI_FILTER$rh100), ONP_CHM_GEDI_FILTER$hoh_2013_dsm_4 , xlab="GEDI rh100", ylab="ALS CHM (ft)",
main="Veg Height")
abline(reg, col="blue")
REQUIRED FIGURES: First screenshot of the filtered GEDI points using mapview and second, a screenshot of the final regression model.
In total for the lab so far you should have:
Now for the last part of the lab
All your figures and tables need to be fully captioned and submitted in one PDF document. All figures need to be high quality and captions need to describe in brief what the image is of, what data was used, and how it was created.
In addition please explain in 300 words or less:
What is lidar?
What is ALS?
What is TLS vs MLS?
What is GEDI?
How has lidar been used for quantifying forest structure and topography?
You will include all of the required figures and tables as well as the short essay, same as 433.
In addition, please explain in 300 words or less: